Revisit on "Ruling out chaos in compact binary systems" 
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Full general relativity requires that chaos indicators should be invariant in various spacetime 
coordinate systems for a given relativistic dynamical problem. On the basis of this point, we 
calculate the invariant Lyapunov exponents (LEs) for one of spinning compact binaries in the 
conservative second post-Newtonian (2PN) Lagrangian formulation without the dissipative effects 
of gravitational radiation, using the two-nearby-orbits method with projection operations and with 
coordinate time as an independent variable. It is found that the actual source leading to zero LEs 
in one paper but to positive LEs in the other does not mainly depend on rescaling, but is due to 
two slightly different treatments of the LEs. It takes much more CPU time to obtain the stabilizing 
limit values as reliable values of LEs for the former than to get the slopes (equal to LEs) of the 
fit lines for the latter. Due to coalescence of some of black holes, the LEs from the former are not 
an adaptive indicator of chaos for comparable mass compact binaries. In this case, the invariant 
fast Lyapunov indicator (FLI) of two nearby orbits, as a very sensitive tool to distinguish chaos 
from order, is worth recommending. As a result, we do again find chaos in the 2PN approximation 
through different ratios of FLIs varying with time. Chaos cannot indeed be ruled out in real binaries. 

PACS numbers: 04.25.Nx, 05.45. Jn, 95.10.Fh, 95.30.Sf 



The chaotic behavior of nonlinear dynamical systems 
has become a very interesting subject in relativistic as- 
trophysics [1] . Especially the dynamics of binary systems 
of spinning compact objects in the frame of general rel- 
ativity does deservedly receive a great deal of attention. 
Merging binaries are regarded as the most promising can- 
didates for future ground- and space-based gravitational 
wave detectors, such as LIGO [2]. The successful detec- 
tion is necessary to rely on the matched filtering tech- 
nique with the theoretical gravitational wave templates 
matched to experimental data containing a lot of instru- 
mental noise. However, chaos in the gravitational wave 
sources would affect the treatment of observational data, 
for example, signals not to be drawn out of the noise. 
For this reason, there have been a series of articles [3-9] 
for discussing whether spinning compact binaries can ex- 
hibit chaos. In the light of the method of fractal basin 
boundaries, an earlier paper [3] emphasized that the con- 
tribution of spin-orbit (SO) and spin-spin (SS) coupling 
is in favor of chaos for the case of comparable mass com- 
pact binaries in the Lagrangian formulation to 2PN or- 
der with the dissipative effects of gravitational radiation 
turned off. While another work [4] suggested ruling out 
chaos by finding no positive LEs of trajectories along the 
fractal basin boundaries of Ref. [3]. At once, as an an- 
swer to this claim, it was reported in Refs. [5,6] that the 
wrong results of Ref. [4] should be owing to the less rig- 
orous calculation of the LEs of two nearby orbits, with 
unapt renormalization time steps adopted. Further some 
orbits with positive LEs were given. 

Obviously, it is very surprise that Ref. [4] and Refs. 
[5,6] employed the same chaos index — the largest LE, 
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but gave completely different dynamics to the same 2PN 
equations of motion for spinning compact binaries. Al- 
though the latter pointed out the problem of the former, 
such interpretation seems still to be ambiguous and puz- 
zling in retrospect. Here are more thorough comments on 
these works. LEs, as a common chaos indicator, measure 
the rate of exponential divergence between neighboring 
trajectories in the phase space. There are two different 
methods to compute them. Historically, the tangent vec- 
tors ^(0) and ^{t) about a given trajectory at times and 
t are used to define the maximum LE: A = limt_>oo x(^)j 
with x{t) = (lA)ln[|4(t)|/||(0)|]. The technique for get- 
ting the LE is called as method 1 (Ml). Usually it is 
a cumbersome task to derive the variational equations 
associated to the tangent vector for complicated dynam- 
ical systems. For an alternative procedure to Ml, a sim- 
pler way, M2, is to adopt the distance d{t) in the phase 
space between two nearby trajectories as an approxima- 
tion to the norm of the tangent vector ^{t) such that 
x{t) = {l/t)\n[d{t)/d{0)]. It is for a suitable choice of 
the starting separation d{0) and of the rescaling interval 
that M2 gives almost the same values of LEs as Ml does 
(for details, see Ref. [10]). As a practical application of 
M2, traditionally one plots a curve of Inx(i) vs Int. A 
negative constant slope of the curve means the regularity 
of the system. If the slope tends gradually to zero and 
In x{t) reaches nearly a stabilizing value, the bounded 
system becomes chaotic. The diagram method is marked 
as M2a. In addition, there is another diagram method 
(M2b) by plotting \n[d{t) / d{0)] vs t. It is vital to per- 
form a least-squares fit on the simulation data to work 
out the slope x of the fit line \n[d{t) / d{0)] = x^, as the 
largest LE. There should have been no difference between 
M2b and M2a in principle, but the former superior to the 
latter lies in that it is much easier to identify the linear 
growth of \u[d{t) / d{0)] than to identify the convergence 
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of Inx(t). In other words, generally it costs a rather long 
time for hix{t) to converge a limit value in the chaotic 
case. As to the fit slope, perhaps it is not very true but 
can always easily be seen even if time is short. In partic- 
ular, the difference is rather explicit when the authors of 
Ref. [4] and those of Refs. [5,6] used M2a and M2b to 
treat the LEs of compact binaries, respectively. Thus we 
think the very slow convergence of Inx(t) leading to the 
"false" LEs in Ref. [4] , but do not agree with the authors 
of Refs. [5,6], who claimed that the wrong LEs in Ref. 
[4] depend on the rescaling. This will be further checked 
in our next numerical experiments, where both M2a and 
M2b adopt the same rescaling interval. 
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FIG. 1: Invariant LEs of the considered three orbits, Fi, r2 
and Fa , by use of the method M3a. 

Now we conclude several problems that appear in Ref. 
[4]. (1) The initial distance problem. An important 
point to note is that relatively large and small initial 
separations are not permitted when M2 is used. For a 
machine double-precision environment with an order of 
10~^^, the starting distance d{0) with a magnitude of 
10~^ is viewed as the best choice [10]. However, Ref. [4] 
used d{0) = 10~^" that must give rise to the overesti- 
mation of LEs if integration time is long enough. The 
rescaling that brings roundoff errors can certainly have 
an effect on LEs, but the initial distance is more impor- 
tant to affect LEs than the rescaling. (2) The integration 
time problem. In general, it is not true for the authors 
of Ref. [4] to declare the absence of chaos in compact 
binaries by finding no stabilizing values of Inx(i) only 
within a time span of a lower limit on the Lyapunov time, 
tx = 1/A, with many times greater than the typical inspi- 
ral time. Perhaps the authors considered that chaos after 
the inspiral time scale will not affect the dynamics of coa- 
lescing compact binaries. It is correct. However, for M2a 
it usually takes many and many times greater than the 
Lyapunov time (rather than the inspiral time) for Inx(t) 
to approach to a certain stabilizing value. For instance. 



Ref. [11] found that x(i) of an orbit in 3-dimensional 
systems seems to have been stabilized to a value near 
0.0005 up to t = 220000 and then abruptly jumps to a 
value around 0.01 up to i = 1600000 (see Fig. 10a in Ref. 
[11]). That is to say, it is completely impossible to arrive 
at the reliable value about 0.01 of LE when the orbit is 
integrated to the Lyapunov time, 100. In sum, the or- 
bits of compact binaries must be integrated numerically 
for sufficiently long times, otherwise there arc unreliable 
LEs. Unfortunately, coalescence does no longer give a 
chance to numerical integration. As an illustration, for 
the conservative system in which gravitational radiation 
is turned off, the coalescence is not a consequence of en- 
ergy loss but just that these chaotic orbits happen to veer 
too close at some stage and merge. It should be possible 
in principle to find pairs that execute enough orbits that 
they do not coalesce before a lengthy integration has been 
performed. (3) The coordinate gauge invariance problem. 
There is a long history of the problem using LEs reliably 
in a curved space. In the mixmaster cosmology there was 
a long standing debate that the LEs were zero therefore 
there was no chaos [12-15]. This was a wrong conclusion 
and was an artifact of the spacetime slicing. Many in- 
dependent groups have been engaged to this field. For 
example, Imponente and Montain [16] gave an invariant 
treatment of LEs by projecting a geodesic deviation vec- 
tor for the Jacobi metric on an orthogonal tetradic basis 
so that they could successfully gain an insight into the dy- 
namics of the mixmaster cosmology. So did Motter [17], 
who addressed directly the issue of the invariance of LEs. 
The invariant LEs in these works arc mainly focused on 
the time evolution of the gravitational field itself. How- 
ever, relativistic compact binaries are attributed to the 
geodesic or nongeodesic motion of particles in a given 
gravitational field. Ref. [3] used fractal basin boundary 
methods to detect chaos in black hole pairs. It should 
be mentioned that the original conclusion that there is 
chaos in spinning binaries was made using a coordinate 
invariant approach. There is no ambiguity in that ap- 
proach. But the fractals can't tell one more details of 
the dynamical features, such as the timescale for chaos 
to set in. Thus, it is fair to try to find a good invari- 
ant version of the LE. There are several points regarding 
this. From the physical point of view, it is questionable 
that the LEs, based on coordinate time t and the Carte- 
sian distance d{t) between two nearby trajectories in the 
phase space of 12 dimensions, are used to discuss the 
dynamics of this system in Ref. [4]. Compact binaries 
are such strong fields that general relativistic effects be- 
come very apparent. On the other hand, general relativ- 
ity admits a free choice of space and time coordinates so 
that the spacetime coordinates usually play book-keeping 
only for events. Therefore, physical observable quanti- 
ties, such as the distance and the time, should be de- 
fined as proper quantities instead of coordinate quanti- 
ties. This is a basic point from the theory of observation 
in general relativity. Following this idea, Ref. [18] used 
proper time r of an "observer" and a proper configura- 
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FIG. 2: Invariant LEs as plots of ln[AL(r)/AL(0)] vs r, based on the method M3b. Curve Fi, with Lyapunov time fx = 
2991M = 10.9To, nearly consists with the lower line of Fig. 4 in Ref. [6]. The integration time is t = 46000M for Fi, while 
t= lO^M for Fa or F3. 
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FIG. 3: The evolution of invariant FLIs with proper time. 



tion space distance AL(t) between the observer and his 
"neighbor" particles to construct an invariant LE (M3): 
A = lim^^oo x(t), where x(t) = (l/r) ln[AL(T)/AL(0)]. 
AL(t) = {haisAx^Ax^y^^, with the space projection 
operator of the observer h"^ = g°'^ + If^U^ , and Ax^ 
being the deviation vector from the observer to the neigh- 
bor. Here g"^ and ?7" stand for metric tensor and 4- 
velocity of the observer, respectively. In practice, M3 is 
no other than a direct modified and refined version of 
M2. Naturally, M2a and M2b are corresponded to M3a 
and M3b. As an illustration, the coordinate time t still 
remains of a common time variable in the equations of 
motion for the two particles, while the proper time r is 



from integration of the equation about dr/dt. For the 
special case of r ~ Inf, M3 fails to work well. In spite 

of that, we do believe that M3 will be very useful and 
simple to investigate spinning compact binaries, because 
it is in the coordinate time t that the equations of motion 
for the compact binaries have been given by Ref. [2] , and 
r and t have no the approximately logarithmic relation 
at all. (4) The power spectra problem. The authors of 
Ref. [4] did not find chaotic behavior in terms of the 
power spectra. This is because the power spectra are dif- 
ficult to distinguish among complicated periodic orbits, 
quasi-periodic orbits and weakly chaotic orbits. Usually 
the power spectra are not recommended to be a criterion 
for evaluating chaos. 

One main aim of the present paper is to re-review the 
results of Ref. [4] , as has been stated above. The other is 
more important to use the covariant chaos indicator M3 
(M3a and M3b) to investigate spinning compact binaries 
so that we take the opportunity to examine the related 
results in some references [4-6]. Considering the slow 
convergence of LEs and the possible coalescence of two 
stars, we suggest adopting a sensitive tool for detecting 
chaos — the invariant FLI of two nearby trajectories in a 
curved spacetime [19]: FLI{t) = logio[AL(T)/AL(0)]. 
The related details and applications of FLIs can be seen 
in Refs. [19-23] . It stretches exponentially with (proper) 
time for the chaotic orbit, but grows linearly with time in 
the regular case. Throughout the work wc use; units c = 
G = 1 and the signature of a metric as (— , -|-, -|-, -|-), and 
let Greek subscripts run from to 3 and Latin indexes 
from 1 to 3. 

In compact binaries, the evolution equations about 
the relative position x and velocity v for body 1 rela- 
tive to body 2 at 2PN order in harmonic coordinates 
are ac = ajv + oipjv + ax.^so + a^pN + a2SS- The 
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numbers and the letters denote the order of the PN 
expansion and type of the contributions to the rela- 
tive acceleration, respectively. The two spins precess 
hy Si = X Si {i = 1,2). Their explicit forms 
can be seen in Ref. [2]. Now let rrij be mass of 
body I, and the total mass M = nii + m2- In addi- 
tion, we specify {y^,Vi) as position and velocity of each 
body in the center-of-mass (CM) frame. The relations 
among three coordinates t/i, 1/2 and x at 2PN order are 
= {m2/M)x + YT_PN{x,v) + YT^,5{Si,S2) + Y2PN{x,v) 
and ^2 = -imi/M)x + Yipn{x,v) + FlsCSi, Sa) + 
Y2pn{x,v) [24]. Here, the IPN and 2PN terms can be 
found in Ref. [25], while the 1.5 order term is given by 
Ref. [26]. On the other hand, the proper time r of body 
1 in the CM frame satisfies the equationdr/di = [—(.900 + 
2goiv\+gijvlv{)]^/'^. gap, as a function of (2/1,1/2; '^i, '"2), 
is the 2PN metric tensor at body 1. Each metric com- 
ponent is made of the related potentials at the location 
of body 1, and each potential is the sum of the non-spin 
piece and of the spin part. The non-spin part is presented 
by Ref. [27] , and the spin piece is listed in Ref. [26] . See 
also Ref. [28] that contains the sum of the two parts. 
As an illustration, the 2.5 order terms in the references 
are dropped. This physically corresponds to dropping 
dissipative terms. 

Clearly, the coordinate time t plays an important role 
in connection with the motion of body 1 and of body 
2, and the relative motion between the two bodies, but 
proper time does not since it differs for each of three mo- 
tions. This gives us a good chance to apply M3 and the 
metric ga0 to study the dynamics of orbits around body 
1 in the CM frame [29] . The implementation is described 
briefly. We integrate the equations (7) of the relative mo- 
tion, the spin equations (8) and the proper time equation 
(10) numerically together by using a fifth-order Runge- 
Kutta-Fehlberg algorithm with an adaptive coordinate 
time step. At once, we can get S'l, S2, x, v, and r, at 
coordinate time t. Then {Vi,V2]Vi,V2) are determined, 
and body 1 has its 4-velocity C/ = (|^,t;J^,u^|^,w?|^). 
Now body 1 is chosen as an observer, who can measure 
the proper distance Ai to his neighboring orbit. Note, 
the neighboring orbit is not the orbit at which body 2 
stays, and is from an orbit nearby body 1, caused by 
a slight separation of the relative position. In a word, 
numerical integration is to carry out in the relative co- 
ordinate system, but the relativistic dynamics is to in- 
vestigate in the CM frame and whether chaos or not is 
measured by body 1. This is entirely different from the 
treatment of other references, where the Newtonian dy- 
namical methods are used to consider the relative motion 
in spinning compact binaries. 

Let us re-calculate the LEs of three orbits that had 
been studied in Ref. [6]. The related initial conditions 

and parameters of the orbits are listed here. Orbit Fi: 
phase space variables {x, v) = (5.5M, 0, 0, 0, 0.4, 0), mass 
ratio P = ma/mi = 1/3, spin magnitudes 5, = m,, and 
spin directions 61 = tt/2 and 62 = 7r/6. Fa: {x,v) = 
(5.0M, 0,0, 0,0.399,0), /3 = 1, 5, = m^, 9^ = 38°, and 



02 = 70°. F3 is the same as T2 but only 0.428 replaces 
0.399. In addition, let only the first component x of the 
initial relative position of each trajectory have a very 
small deviation. Ax = 10~^M, then we get its corre- 
sponding neighboring orbit. Following M3a, we draw 
plots of logio xiT) vs logio about the LEs of the three 
trajectories in Fig. 1. They all drop before proper time 
r spans lO^M « 3636To (T^ = 275M, the average period 
of orbit Fi). For Fi, the LE time looks to get a reliable 
value, T\ = 4675M = 17.0To. It is larger than the value 
tx = 3080M = II.2T0 given in Ref. [6]. Perhaps one ad- 
dresses a question whether x can still remain the value 
of 1/t\ if numerical integration continues. We have no 
way to answer it since numerical integration has to end 
after t = 1.37 x lO^'M (or r = 1335980M), when the 
two objects coalesce. Above all, neither Fa nor has 
any acceptable stabilizing value when integration time 
t reaches lO^M. Additionally, we did not find any dif- 
ference from these results by making several tests with 
different rcnormalization time steps. This seems to show 
that the results in Ref. [4] are reasonable. However, the 
case is completely different when M3b is used. Seen from 
the calculations including the rescaling, M3b is almost 
the same as M3a, but only a slight difference between 
them lies in plotting ln[AL(r)/AL(0)] vs t instead of 
plotting logio x{t) vs log^g r. Another point to note is 
that M2b {rather than M3h) without rescaling was used 
in Ref. [7], but our M3b employs rescaling. In addi- 
tion, there is a difference that the authors of Ref. [7] use 
the Hamiltonian formulation in ADM coordinates and 
not the Lagrangian formulation in harmonic coordinates. 
Although the two approaches are approximately related, 
they are not exactly equal. For instance, the constants of 
motion are exactly conserved in the Hamiltonian formu- 
lation, while they are approximately in the Lagrangian 
formulation. Ref. [7] also works to 3PN order. As shown 
in Fig. 2, it takes no long enough time to explicitly see 
the presence of positive slopes of the fit lines for Fi and 
Fa, but to do that of about zero slope of the fit line for 
F3. This means chaos of Fi and F2, while the regular- 
ity of F3. It is what can be seen in Refs. [5,6]. It is 
sufficiently argued that the LEs converge much faster for 
M3b than for M3a. As an illustration, t\ is more reliable 
than T\. This is because the longer numerical integra- 
tion becomes, the more accurate the values of LEs are. 
In fact, the fit slope (its inverse being 4680M) of Fi in 
Fig. 2 is very close to the LE of Fi in Fig. 1 when 
integration times are the same. Now, we can say quite 
plainly that a reliable conclusion is that there is chaos 
in the 2PN system. So can the authors of Ref. [7] , who 
have already confirmed the existence of chaos in the 2PN 
Hamiltonian formulation through positive LEs. As men- 
tioned above, although the LEs converge much faster for 
M3b than for M3a, long integration times are still needed 
to get reliable values of LEs even if M3b is considered. 
Noting this, we recommend to use a quicker indicator, 
the invariant FLI given by Eq. (6). Its algorithm can be 
found in Ref. [19]. Fig. 3 displays that FLIs of Fi and F2 
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increase exponentially with logj^Q t, but that of does 
algebraically. Thus Fi and r2 arc chaotic, (chaos of Fi is 
much stronger than that of F2) but F3 becomes ordered. 
It is worth emphasizing that the three orbits can be dis- 
tinguished clearly in practice when proper time adds up 
to lO^M. Consequently, the onset of chaos in the 2PN 
Lagrangian approximation is proved again through dif- 
ferent ratios of FLIs varying with time. 

The summary is included as follows. For conceptual 
clarity, it is necessary to apply chaos indicators indepen- 
dent of the choice of coordinate gauge to analyze the 
dynamics of relativistic gravitational systems. Since co- 
ordinate time is a good medium in connect with the mass 
centered motions of both body 1 and body 2, and the rel- 
ative motion in spinning compact binaries, we think that 
M3, as an invariant indicator, is a good tool to study 
these systems. Using M3, we estimate the LEs on the 
mass centered motion of body 1 rather than on the rela- 
tive motion considered by other references. We find that 
the orbits must be calculated for long enough times in 
order to get stabilizing limit values as reliable LEs for 
the case of comparable mass compact binaries. On the 
other hand, we track that the exact source of both the 
failure of Ref. [4] and the success of Refs. [5,6] in the 
computation of LEs docs not stem from the rescaling, 
but is based on two slightly different treatments of LEs, 
M3a and M3b. At most cases, the LEs converge much 
faster for M3b than for M3a. However, coalescence of 



the black holes makes it impossible in some cases to have 
enough numerical integration. This shows that M3a is no 
longer a suitable indicator to quantify chaos in spinning 
compact binaries. Additionally, it should be noted that 
although it is rather easier to get the LEs for M3b than 
for M3a, long integration times are still needed to get re- 
liable values of LEs when M3b is adopted. In this sense, 
the invariant FLI in a curved spacetime is a very fast and 
valid technique to detect chaos from order. Still, a reli- 
able conclusion is that there is chaos in the conservative 
2PN Lagrangian system. Of course, the 2PN approxima- 
tion is so poor that there has been left an open question 
whether real binary systems with better approximations 
exhibit chaos [30] . Saying this another way, now one does 
not say that chaos can be ruled out in real binaries. In 
future, we will discuss a wider application of the FLI in 
detailedly investigating the dynamics of spinning com- 
pact binaries. 
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